Bright Matter- Wave Soliton Collisions in a Harmonic Trap: Regular and Chaotic Dynamics 
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Collisions between bright solitary waves in the ID Gross-Pitaevskii equation with a harmonic potential, which 
models a trapped atomic Bose-Einstein condensate, are investigated theoretically. A particle analogy for the 
solitary waves is formulated, which is shown to be integrable for a two-particle system. The extension to 
three particles is shown to support chaotic regimes. Good agreement is found between the particle model and 
simulations of the full wave dynamics, suggesting that the dynamics can be described in terms of solitons both 
in regular and chaotic regimes, thus presenting a paradigm for chaos in wave-mechanics. 
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The presence of chaos in quantum systems is a topic of in- 
tense theoretical and experimental interest fl]. A signature 
of classical chaos is the ergodic filling of regions in phase 
space. Applying this criterion in the search for chaos in wave- 
mechanical systems, for example the linear Schrodinger equa- 
tion in quantum mechanics, one confronts the uncertainty re- 
lations, which dictate that trajectories are smeared out. In par- 
ticular, chaos is impossible to observe when dispersion domi- 
nates over the exponential divergence of neighboring trajec- 
tories. For this reason, non-dispersive waves such as soli- 
tary waves or solitons are of particular interest in the study of 
chaotic dynamics. In this case, particle-like chaotic behaviour 
may be well-defined in wave-mechanical systems. 

Solitary waves may be found in solutions to nonlinear wave 
equations where the nonlinearity counteracts the dispersion of 
a wave-packet such that it retains its form as it propagates. 
The nonlinear Schrodinger equation (NLSE) is such an equa- 
tion, employed to describe diverse physical systems; notably 
light propagating in fibre-optics ||2|], and more recently as an 
approximation to the many-body dynamics of dilute atomic 
Bose-Einstein condensates (BEC) |3], where it is called the 
Gross-Pitaevskii equation (GPE). Solitons are solitary waves 
that emerge unscathed from collisions, up to shifts in posi- 
tion and phase |4]; this behaviour is reminiscent of particle 
behaviour, motivating the particle-like name soliton. The ho- 
mogeneous ID NLSE with attractive nonlinearity supports 
bright soliton solutions |4], so-called because they represent 
a peak (rather than a trough) in intensity in a nonlinear-optical 
setting, or in particle density in BEC. Macroscopic quantum 
states consisting of multiple bright matter- wave solitary waves 
present an interesting testing ground for wave chaos. In gen- 
eral all wave-packet evolution predicted by the Schrodinger 
equation is periodic or quasiperiodic, due to its linearity. The 
nonlinearity in the ID GPE and associated solitary wave so- 
lutions may break all periodicity, leading to the realization of 
ergodic behaviour, i.e., recognizably chaotic dynamics. 

Bright solitary waves have been the subject of substantial 
experimental and theoretical investigation in the context of 
nonlinear optics QllBl, and BEC O S S 113 U?] ■ Notably, 
chaotic and regular soliton behaviour have been observed the- 
oretically in a NLSE with a 5-kicked rotor potential 1 12]. In 
BEC experiments, the magnetic or optical trap employed to 



confine the constituent atoms introduces a position-dependent 
potential. In this Letter, we investigate to what extent solitary 
wave collisions in a harmonic potential provide a paradigm 
for particle-like chaotic behaviour in a wave-mechanical sys- 
tem. Due to this potential, bright solitary waves in a harmon- 
ically trapped system are not true solitons; however, it will be 
shown that the particle nature of the solitary waves is very pro- 
nounced, so the bright solitary waves in this system will from 
now on be called solitons. In order to test the extent of the soli- 
ton behaviour, we introduce a particle model, adapted from a 
nonlinear optics context for a NLSE with a sinusoidal external 
potential |5]. In this model, constructed for the regime where 
the solitary waves are well-separated before and after colli- 
sions, the waves are modeled as interacting classical particles. 
Within this model, we show that the two soliton case is inte- 
grable, but for three (or more) solitons one can expect chaotic 
dynamics. The results are compared to numerical solutions of 
the NLSE, and provide a probe of the coexisting particle and 
wave properties of bright solitary waves. The most surprising 
result is that the particle-like behaviour is preserved even in 
the chaotic regime. In contrast to the linear Schrodinger equa- 
tion, where the evolution of localized wave-packets is rapidly 
disrupted in regimes supporting classically chaotic dynamics 
1 1], the soliton solutions appear to be remarkably robust. At 
temperatures encountered in atomic BEC experiments (nK or 
lower), the atom-atom interaction potential may generally be 
replaced by an effective contact interaction, quantified by the 
i-wave scattering length, a. Depending on species, this may 
be positive or negative. By exploiting a Feshbach resonance, it 
may also be tuned using an external magnetic field 1 13]. In the 
case of a trapped, almost fuUy Bose-condensed dilute atomic 
gas, the dynamics are largely governed by the following GPE: 
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(1) 
where A^ is the total number of atoms, m the atomic mass, 
and ^(r, t) the condensate mode-function, normalised to one. 
The atom-atom interactions are quantified by gT,D - A-ntfa/m, 
where, in this Letter, a is negative. The proportion of non- 
condensate atoms is thus assumed to be negligible; it can be 
shown, however, that linear instabilities in the GPE directly 
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FIG. 1: (Colour online). Trajectories in the particle model (yellow 
lines) plotted over density distributions predicted by ID dynamics in 
the homogeneous GPE. The trajectories correspond to solitons col- 
liding with a relative phase of the golden ratio = (1 + V5)/2 in (a) 
and (c), and a relative phase of = ;r(l + V5)/2 in (b) and (d). In (a) 
and (b) the incoming speeds of the solitons are -0.1 |gio|A'/S and 0.2 
IgiolA'/S, and in (c) and (d) the incoming speeds are -0.1 IgiolA'/ft 
and 0.3 IgiolA'/^. Taking the parameters of the system to agree with 
recent experiment |8], the solitons have equal effective masses, the 
axial trapping frequency is 10 Hz, and other parameters (radial trap 
frequency of 800 Hz, atomic mass and scattering length of 'Li, and 
5000 particles per soliton). The unit of x is then equal to 3.6 yum, and 
a unit of f to 1.4 ms. 



imply Ill4tll5r that the population of the non-condensate com- 
ponent may rapidly become significant. Regimes where soli- 
ton collision dynamics have a chaotic character are thus of ad- 
ditional interest, as they may coincide with a greater tendency 
for linear instability, and hence implicitly with condensate de- 
pletion 1 16]. 

The regime of interest for the study of solitons is the quasi- 
ID case, where the atoms are trapped in a radially tight har- 
monic trap with loose harmonic axial confinement. We may 
assume a harmonic ground state (Gaussian) ansatz in the ra- 
dial direction, since the harmonic potential energy dominates 
over the interaction energy. The GPE then takes the form: 
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where x is now measured in units of fi^/mlgiolA^ and t in units 
of Tv' jm\g\i)^N'^ with giu - IfiUra, Ur is the radial trapping 
frequency and u is the axial frequency in our units of inverse 
time. In the case of vanishing harmonic potential, an exact 
solution exists |4], comprising an arbitrary number of well 
separated solitons taking the form: 
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where qi - vit + xot is the position of the peak of the /th soli- 
ton; XQi is the peak position at f = 0; ao/ - vixoi is the phase 
at X = 0, f = 0; and v, are the soliton velocities. Our nor- 
malisation condition implies that 2 ■ ' 4?7,- - 1, where A^^ is the 
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FIG. 2: (Colour online). Poincare section of the three-soliton system 
with H=10. Regions corresponding to trajectories in figures 3(a) to 
3(c) are labeled, and highlighted using larger, darker points. The sec- 
tion corresponds to the momentum p,- and position q^ of the "asym- 
metric stretch" mode when the "stretch" mode coordinates q^ = 0, 
Pc < 0. The figures correspond to the regime where the solitons have 
equal effective masses, the axial trapping frequency is 10 Hz, and the 
other parameters (radial trap frequency of 800 Hz, atomic species 
mass and scattering length of 'Li, and 5000 particles per soliton) 
correspond to recent experiment Ql. 



number of solitons present. When these solitons emerge from 
collisions, they suffer position shifts dependent on their initial 
speeds v,, and the quantities //, known as effective masses only. 
The outgoing soliton motion does not depend on the relative 
phase |4]. This property is illustrated in Fig. 1 where we com- 
pare the solution of the ID GPE with a particle model which 
is outlined below. Note that the density distribution during the 
collision does depend on the relative phase, but this does not 
influence the outgoing trajectories. 

The particle model follows the approach of Scharf and 
Bishop 1 5], which reproduces the position shifts following 
collisions, and also reproduces the motion due to the trap- 
ping potential, whilst neglecting the phase behaviour. This ap- 
proach is appropriate for regimes where the solitons are well 
separated between collisions. This approach is not appropri- 
ate for soliton trains, as observed by Strecker et al. |8|], and 
modeled by Gordon [IT] in the absence of any external po- 
tential, and by Gerdjikov et al. 1 181 in the case of harmonic 
confinement, where the solitons are never well separated, and 
the phase difference has an important effect. Parker et al. have 
modelled bright matter-wave soliton collisions 1 1 1 ] using the 
3D GPE at or near the quasi- ID regime. Taking the full 3D 
GPE dynamics into account highlights some important devi- 
ations from the ID dynamics: in particular, collapse may oc- 
cur during collisions of soUtons having slow approach speeds, 
with sensitivity to the relative phase of the solitons. However, 
above a particular threshold velocity, the quasi- ID model con- 
sidered in this Letter can be expected to hold, including the 
observed phase-insenstivity of the soliton collision dynamics. 

Following reference |5], the action of the external potential 



is deduced by using the one soliton solution of the homoge- 
neous case as an ansatz for the system with a harmonic poten- 
tial, and evaluating the norm and energy functionals, which 
are known to be constant. This gives equations of motion for 
the solitons in the external potential. Added to this is an inter- 
particle potential which reproduces the position shifts of the 
solitons on emerging from collisions with each other, inferred 
from the exact solution to the homogeneous NLSE, which 
are assumed not to change upon the addition of the position- 
dependent external potential. The Hamiltonian for an arbi- 
trary number of solitons (A^j) is given in the particle analogy 
by: 



the reduced system Hamiltonian 






'1 277, 



- ^ 277,77/ 77; + 77^)sech^ 



1 <i<j<Ns 



rii + rjj 



iqi - qj) 



(4) 



in which the reason for the interpretation of 77,- as effective 
masses is clear. This Hamiltonian models the positional dy- 
namics of the soliton peaks. In the case of two solitons (A^, - 
2) with identical effective masses, the following coordinates 
may be defined: the centre-of-mass position Q :- (qi + ^2)/2 
and the relative position q :- qi - qi- The Hamiltonian [Eq. 
01 reduces to 
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where P is the momentum conjugate to Q, and p the momen- 
tum conjugate to q. The Hamiltonian is seperable into two 
parts: the centre-of-mass energy E (dependent on P and Q 
only), and the "relative energy" e (dependent on p and q only). 
There are thus two independent constants of the motion, E and 
e, as many as there are degrees of freedom. Hence, the particle 
model for two solitons is integrable and the dynamics must be 
completely regular The same argument holds in the simple 
generalisation to non-identical masses. 

In the case of three solitons (A^, = 3), the situation is dif- 
ferent. When the masses are identical, a useful coordinate 
change may be made to QtIt] '■- {q\ + q2 + ?3)/3, the centre- 
of-mass position, and qc/j] := (^1 - q3)/2 (corresponding 
to the "stretch" mode) and ^,777 :- (qi + q^ - 2q2) (corre- 
sponding to the "asymmetric stretch" mode), the normal co- 
ordinates of the system for small displacements from the ori- 
gin. The stretch modes are similar to the vibrational modes 
in a tri-atomic molecule 119]; as the system is constrained to 
ID, however, there is no analogue of the molecular bending 
mode. Rescaling time I - jft , and introducing the momenta 
Pc = 2qc and p,- = q,/€>, we may remove the centre-of-mass 
behaviour from the problem, as it decouples from the other 
degrees of freedom. (An equivalent treatment is possible for 
non-identical masses, 7/,-.) The resultant Hamiltonian we call 
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This Hamiltonian, describing the two remaining degrees of 
freedom, is not separable, and it is necessary to integrate the 
equations of motion numerically to analyse the system's be- 
haviour. Poincare sections illustrate regions of regular and 
chaotic behaviour Figure |2 shows a Poincare section corre- 
sponding to the momentum pr and position q^ of the "asym- 
metric stretch" mode when the "stretch" mode coordinates 
qc - 0, Pc < 0; other sections would be equally illustrative 
of the qualitative behaviour. The behaviour is regular at large 
positive values of H, but as H is reduced, chaotic behaviour 
emerges, characterised by ergodic regions in between regular 
tori. For small H the system is mostly an ergodic sea, with is- 
lands of stability. A regime with both ergodic regions and reg- 
ular tori is plotted in Fig. |2 Consideration of the form of the 
reduced-system Hamiltonian [Eq. (|6jl] shows that without the 
interaction the system is integrable, as it becomes a decoupled 
pair of harmonic oscillators. When H is large and positive, the 
interaction part of the Hamiltonian (which is always negative) 
should give a relatively small contribution to the Hamiltonian, 
compared to the integrable part of the Hamiltonian (which is 
always positive). When H is reduced, this is no longer the 
case, and chaos is emergent. 

Figure |3 shows a comparison of trajectories in the parti- 
cle model with results from integrations of the ID GPE [Eq. 
0]. The comparisons illustrate the good agreement between 
the particle model and the ID GPE in the regimes in which 
the particle model is valid, i.e., when solitons are well sep- 
arated between collisions [Figs. |3a) and|3b)], even when 
the motion is chaotic [Fig. |3jb)]. When two of the solitons 
are not well separated [Fig. |3jc)], the ID GPE simulation 
shows that a "bound state" is formed, which looks like a single 
"higher-order" soliton with an excited breathing mode. The 
particle model does not predict well the behaviour within the 
"bound state", but does give a good prediction of the centre- 
of-mass motion of the "bound state" and its interactions with 
the other soliton; it is likely that the behaviour of the density 
of the "bound state" is strongly coupled to the phase behaviour 
within the "bound state". 

These results show that harmonically trapped solitons have 
strong particle characteristics. Even in chaotic regimes, where 
the exponential growth of linear instabilities is most preva- 
lent, the soliton solution is remarkably robust. An echo of 
the wave-equation origin of the particle model is that, consis- 
tent with the attractive interaction potential, the particles pass 
through each other subsequent to collisions. The payoff of us- 
ing the particle model is that it is quicker to solve than the ID 
GPE (involving four real variables, rather than a continuous 
complex- valued function); thus Poincare sections (see Fig.|2li 
can rapidly build up a quaUtative idea of the many-soliton be- 
haviour 
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FIG. 3: (colour online). Trajectories in the particle model (yellow 
lines) plotted over density distributions predicted by ID GPE dynam- 
ics, corresponding to (a) a regular orbit, (b) a chaotic orbit, and (c) a 
bound state, as shown in Fig.|2{b). In the centre of the bound state, 
the density increases to approximately 0.5 normalised units, but our 
scale is pinned at 0.14 units in order to resolve the low density re- 
gions better. The parameters of the system are H=10, the solitons 
have equal effective masses, the axial trapping frequency is 10 Hz, 
and other parameters (radial trap frequency of 800 Hz, atomic mass 
and scattering length of 'Li, and 5000 particles per soliton) corre- 
spond to the recent experiment Ql. The unit of x is then equal to 2.4 
fira, and a unit of / to 0.6 ms. 



The experimental demonstration of such chaotic dynamics 
in a wave mechanical system requires a relatively straightfor- 
ward adaption of recent experiments on bright matter wave 
solitons 1 8, 9, 10]. For example, a system of 3 solitons can be 
created reproducibly by careful choice of the initial conditions 
llOll . Manipulation of the optical trapping potential during the 
creation of the solitons will allow the solitons' initial veloci- 
ties to be chosen. The chaotic regions of phase space may be 



probed by measuring the sensitivity of the subsequent evolu- 
tion of the density distribution to the initial condition. 

In conclusion, an effective classical particle model has been 
derived for many solitons in the NLSE with a harmonic po- 
tential. This applies to a dilute BEC of attractive atoms in 
the quasi-lD limit of a cylindrically-symmetric cigar-shaped 
trap. Within this model two-soliton dynamics are fully in- 
tegrable and regular, but three solitons may display chaotic 
dynamics when atom-atom interactions are significant. The 
particle model exhibits good agreement with the ID GPE in 
the regime of large separation of the solitons before and after 
collisions, even when the particle motion is chaotic. This con- 
firms the surprising robustness of bright matter- wave solitons, 
as has also been observed experimentally fo, 9, lOV There 
is a good degree of agreement even when "bound states" are 
modeled (states not in a regime of large separation). Chaotic 
regions may also be a useful predictor of regimes of conden- 
sate instability, which can be explored with a fuller treatment 
of the condensate and non-condensate atoms llSll. 



We thank J. Brand, S. L. Cornish, T. S. Monteiro and N. G. 
Parker for useful discussions, and acknowledge support from 
the UK EPSRC. 



[1] M. C. Gutzwiller, Chaos in Classical and Quantum Mechanics 
(Springer- Verlag, New York, 1990); L. E. Reichl, The Tran- 
sition to Chaos In Conservative Classical Systems: Quantum 
Manifestations (Springer- Verlag, New York, 1992). 

[2] H. A. Haus and W. S. Wong, Rev. Mod. Phys. 68, 423 (1996). 

[3] F. Dalfovo etal. Rev. Mod. Phys. 71, 463 (1999); C. J. Pethick 
and H. Smith, Bose-Einstein Condensation in Dilute Gases 
(Cambridge University Press, Cambridge, 2002); L. Pitaevskii 
and S. Stringari, Bose-Einstein Condensation (Clarendon Press, 
Oxford, 2003). 

[4] V. E. Zakharov and A. B. Shabat, Zh. Eksp. Teor. Fiz. 61, 
118 (1971); Sov. Phys. JETP 34, 62 (1972); L. D. Faddeev and 
L. A. Takhtajan, Hamiltonian Methods in the Theory of Solitons 
(Springer- Verlag, Berlin Heidelberg, 1987). 

[5] R. Scharf and A. R. Bishop, Phys. Rev. A 46, R2973 (1992); 
Phys. Rev. E 47, 1375(1993). 

[6] R. Scharf, Chaos, Soliton. Fract. 5, 2527 (1995); H.-H. Chen 
and C.-S. Liu, Phys. Rev. Lett. 37, 693 (1976); N.-C. Panoiu 
etal. Chaos 10, 625 (2000). 

[7] L. D. Carr and J. Brand, Phys. Rev. Lett. 92, 040401 (2004); 
U. Al Khawaja et al, ibid. 89, 200404 (2002); L. Salasnich, 
A. Parola, and L. Reatto, Phys. Rev. A 66, 043603 (2002); 
V. Y. F. Leung, A. G. Truscott, and K. G. H. Baldwin, ibid. 66, 
061602 (2002); V. M. Perez-Garcia, H. Michinel, and H. Her- 
rero, ibid. 57, 3837 (1998); L. D. Carr, C. W. Clark, and W. P 
Reinhardt, ibid. 62, 063611 (2000). P V. Elyutin e/ a/. , Phys. 
Rev. E 64, 016607 (2001); Y S. Kivshar, T. J. Alexander, and 
S. K. Turitsyn, Phys. Lett. A 278, 225 (2001); 

[8] K. E. Strecker et al. Nature 417, 150 (2002). 

[9] L. Khaykovich etal. Science 296, 1290 (2002). 
[10] S. L. Cornish, S. T. Thompson, and C. E. Wieman, Phys. Rev. 

Lett. 96, 170401 (2006). 
[11] N. G. Parker et al, cond-mat/0603059 (unpublished). 
[12] R Benvenuto et al, Phys. Rev. A 44, R3423 (1991). 



[13] S. Inouye et ai. Nature 392, 151 (1998); T. Kohler, K. Goral, et al, Phys. Rev. A 62, 023612 (2000); Y. Castin and R. Dum, 

and R S. Julienne, cond-mat/0601420 (unpublished). Phys. Rev. Lett. 79, 3553 (1997). 

[14] N. N. Bogoliubov, J. Phys. (USSR) 11, 23 (1947). [17] J. R Gordon, Opt. Lett. 8, 596 (1983). 

[15] Y. Castin and R. Dum, Phys. Rev. A 57, 3008 (1998); C. W. [18] V. S. Gerdjikov et al, Phys. Rev. E 73, 046606 (2006). 

Gardiner, ibid 56, 1414 (1997). [19] J. Tennyson and S. C. Farantos, Chem. Phys. 93, 237 (1985). 

[16] S. A. Gardiner, J. Mod. Opt. 49, 1971 (2002); S. A. Gardiner 



